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Abstract. We investigated hydrodynamical post-Newtonian models of selfgravitating 
stationary black-hole-disk systems. The post-Newtonian scheme presented here and also in 
our recent paper is a continuation of previous, purely Newtonian studies of selfgravitating 
hydrodynamical disks rotating according to the Keplerian rotation law. The post-Newtonian 
relativistic corrections are significant even at the 1PN level. The 1PN correction to the angular 
velocity can be of the order of 10% of its Newtonian value. It can be expressed as a combination 
of geometric and hydrodynamical terms. Moreover, in contrast to the Newtonian Poincare- 
Wavre theorem, it depends both on the distance from the rotation axis and the distance from 
the equatorial plane. 

In the technical part of this note we derive virial relations valid up to 1PN order. We show 
that they are indeed satisfied by our numerical solutions. 


1. Introduction 

In a recent paper [Ij we investigated post-Newtonian models of selfgravitating gaseous disks that 
rotate according to the Keplerian rotation law. The analysis presented there is a continuation 
of our previous studies, where such disk systems were investigated in Newtonian theory [ZU31- 
In the Newtonian framework we posed the following problem: Suppose one observes a 
selfgravitating stationary gaseous disk around a central object (modeled by a point-mass) that 
rotates according to the Keplerian rotation law, that is with the angular velocity Uq = u;o/r 3//2 , 
where r is the distance from the rotation axis, and ujq is a constant. For the disk consisting of 
test particles we have uj o = \/GM c , where M c is the mass of the central object, and G is the 
gravitational constant. What does the observed value of ujq correspond to in the case where 
the mass of the disk is comparable with M c ? Is it still the central mass M c , the sum of the 
two masses, or some nontrivial combination of them? It turns out that the selfgravity speeds 
up the rotation of the disk—it rotates faster than this would follow from the Keplerian formula 
involving the central mass M c only. Moreover, the way in which ujq depends on the central 
mass and the mass of the disk is prescribed by the geometry of the disk. Thus, in principle, it 
is possible to measure the masses of Keplerian disks, whenever their geometry is known. This 



procedure was applied to the accretion disk in the AGN of NGC 4258, where the Keplerian 
rotation curve was measured in the maser emission [3]. 

In [T] we extended the Newtonian analysis to the first post-Newtonian approximation (1PN). 
Selfgravitating stationary gaseous disks were investigated before in full relativity (cf. [31 [5]). We 
decided on the post-Newtonian scheme, because of its conceptual simplicity. In particular, the 
notion of the Keplerian rotation has a clear meaning in the post-Newtonian scheme. 

The main result obtained in the 1PN approximation is that the angular velocity profile is 
affected in two different ways—some parts of a disk can be speeded up and the others slowed 
down. Furthermore, the sum of the Newtonian and post-Newtonian components of the angular 
velocity is not anymore a function of the cylindrical radius only, but in general it depends on 
radial and vertical coordinates [I]. 

In this paper we sketch briefly the main equations that constitute the 1PN model and then 
derive virial-type relation that can be used to test the obtained numerical solutions. Suitable 
virial tests valid in the Newtonian case were presented in [6] and [2j. We discuss them here for 
clarity. The post-Newtonian virial identities given here are new and have not been discussed 
in [lj. In the last section of this paper we also show that they are satisfied by our numerical 
models with the accuracy similar to that of Newtonian solutions. 


2. Description of the model 

Our 1PN black-hole-disk models are constructed assuming the metric of the form 
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where we use Cartesian coordinates x° = ct,x 3 = x,x 2 = y,x 3 = z, and c is the speed of light. 

We write the energy-momentum tensor as 

T a P = Mc f ,u BH u B H (5 ( x j + ^ c 2 + h)u a U P+pg a P, (2) 

y /9 u bh 

where the first component describes the point particle (it is proportional to the Dirac delta 
distribution and models the central black hole) at rest, located at the origin of the coordinate 
system; the second one is the energy-momentum tensor of the disk matter. Here M c denotes the 
mass of the point particle; g is the determinant of the metric g = — det(g^ w ). The four-vectors 
and u a denote the four-velocities of the central point-mass and the fluid, respectively. The 
symbol p denotes the baryonic rest-nrass density, h is the specific enthalpy, and p is the pressure. 
In the following sections we will also work with the three-velocity, defined as v l = cu l /u°, 
i = 1,2,3. 

In the remaining part of the article we use standard cylindrical coordinates ( r,z,q i>). We 
consider a stationary, selfgravitating, axially and equatorially symmetric polytropic disk, 
rotating around a central point mass M c according to the Keplerian rotational law = ixor~ 3 ^ 2 . 
We assume that the disk is geometrically bounded by the inner and outer radius r- in and r out , 
respectively. We introduce the notation according to which any quantity £ (if it is necessary) 
is separated into its Newtonian £o and post-Newtonian £i part according to the general pattern 
£ = £o + £i/c 2 - Following [7] we derive basic equations which, separated into their Newtonian 







and post-Newtonian parts, read: 


(3) 

(4) 


A U 0 = AirG (p 0 + M c S(x )), 

ho = -U 0 + J dr(v$) 2 r+ C 0 , 

AA^ = 2 ——- — WirGr 2 poVq , ( 5 ) 

AC/r = 4^G(M c C/ 0 D (0)5(x)+ Pl + 2p 0 + /9o(^o-2C/o + 2r 2 (4) 2 )), (6) 

hi = -C/i - vtAf + 2h 0 (v$) 2 r 2 - J dr(v*) 4 r 3 - h 2 - Ah 0 U 0 - 2U 2 - C u (7) 

where A denotes the flat Laplacian with respect to coordinates (a; 1 , a; 2 , a; 3 ), Co and C i are 
integration constants, and Uq* is the gravitational potential due to the disk only, i.e., Uq = 
—GM C / |x| + U'q . They follow from the conservation law, V a T Q ^ = 0, the continuity equation 
V Q ( pu a ) = 0, and Einstein equations 


Rfiu ndfj.i' — 37 t ^( 8 ) 
A C 

where R^ u is the Ricci tensor and R denotes the Ricci scalar. 

The above system of equations is closed by assuming an equation of state. Our numerical 
solutions are obtained for a polytropic equation of state of the form p = Kp 1 , or equivalently 
ho = K'y/ipy — l)pg 1 (for the 0 th order solution) and h\ = (7 — l)hoPi/po (for the 1PN part), 
where K and 7 > 1 are constants. 

Equations mm should be solved with respect to the Newtonian gravitational potential 
Uo(r,z), the post-Newtonian gravitational potential Ui(r,z), the rotational potential A^r, z) 
and the Newtonian and post-Newtonian enthalpy ho(r,z) and h\(r,z). Any other quantity (the 
density p(r,z), the pressure p(r,z), etc.) can be obtained from these five functions. 

Numerical solutions are obtained as follows. We use the classic Self-Consistent Field (SCF) 
scheme (cf. [8]) to solve the set of Eqs. (|3|) and Eq. ([3]). Given the Newtonian potential Yq and 
the enthalpy ho, we can solve Eq. ((SJ) at once. Finally, we use again the SCF scheme to solve 
the set of Eqs. ([6]) and ([7]). For the Keplerian rotation law the above method converges for all 
values of the parameters that we have tested. 


3. Results 

The main result of our post-Newtonian scheme is the correction to the Newtonian angular 
velocity profile. The well-known theorem by Poincare and Wavre states that Newtonian 
stationary barotropic disks (or stars) rotate with the angular velocity that can depend on the 
distance r from the rotation axis only [9]. It turns out that already the 1PN correction can be 
significant and in general it depends also on z. It can be expressed as 

vf (r, z) = —— t jd r v$ + 2rh 0 d r v$. (9) 
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Note that the above formula involves both geometric and hydrodynamical factors. 

Acceptable 1PN solutions should satisfy the following (quite stringent) conditions: i) 
1 3> \Uo\/c 2 |Gi|/c 4 , ii) 2GM c /c 2 <C ri n , in) c ^$> c s , where c s is a speed of sound. Sample 

numerical models that do satisfy the above conditions are presented in [lj. It turns out that the 
1PN correction to the velocity can be of the order of 10% of the Newtonian value. The reader 
interested in the details of these model may consult pQ. In the remainder of this paper we focus 
on the construction of virial tests that can be applied to our numerical scheme. 




4. Virial identities 

In this section we use Cartesian and cylindrical coordinates. It is implicitly assumed that Latin 
indices refer to Cartesian coordinates. 

The virial identity that can be used to test the Newtonian solution (including the central 
point-mass) was obtained in [ 6 j. It reads 


( 10 ) 


£ P ot T 2-EiQ n + 2_E4h erm — 0, 


where E pot = \ f R3 d?x p(Uq — GM C / |x|) is the total potential energy, Skin = \ / R 3 d 3 x pvQiV l 0 is 
the bulk kinetic energy and S t h er m = I / R 3 d 3 x p is the internal thermal energy. We assume, as 
a virial test parameter, the value 


E pot T 2Ski n T 2-Etherm 

-Spot 


In order to obtain the post-Newtonian relations we rewrite Eqs. ([5]) and <[ 6 |) in a slightly 


different form. Equation Eqs. Q can be written in Cartesian coordinates as 


A A 1 = -16irGp 0 v l 0 . 


Equation © is split in two parts: 


A U[ 

A U'{ = 4ttGM c U® (0)5(x) 
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( 12 ) 


where U\ = U[ + U". The solution for U” is 



GM c L 0 d ( 0) 
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Consider a vector 



Its divergence reads 


did' 



(13) 


For a finite disk (po of compact support), Af. tends to zero sufficiently fast, and 


|x| 2 a* —>• 0, as |x| = \J x 1 + y 2 + z 2 —>• 00 . 

Thus, by integrating Eq. C3D over R 3 , and making use of the Gauss theorem, we see that 



Integrating by parts, one can get rid of the term with diA &. This yields 







where only 0 th order terms are differentiated. 

The above relation can be also written in terms of the vector components in cylindrical 
coordinates. Because of symmetry assumptions, we have A k v q = A^Vq, and 

x l d t (po^o) A k = x l Vi (po't’o) A k = PqA^Vq + x l di (p 0 t$) A <t>- 

This yields 

0 = J ^ d 3 x (^—^pqA^Vq - x l di (povfj A <t> 

The virial relation following from Eq. m can be obtained in a similar way. It is enough to 
consider the divergence of the vector 

) d l U[ - ^diUidfUi 


b l = x l diU[ + -u[ 



It yields 


dib 1 = x l d L u[ + -u[ A u[, 


and analogously 

0 = d 3 x (x l d,U[ + (p, + 2 p 0 + po (ho - 2U 0 + 2 r 2 (^) 2 )) . 


Many different forms of the above relation can be obtained by ‘playing’ with Eq. (ED- a 
helpful relation that can be used here is 


-x'diU? = U'{. 


In the analogy to the Newtonian case we choose as virial test parameters e' = \{e ! 0 + e' b )/e' a \ 
and e " = |(e" + e ")A"|, where 

4 = / d 3 x Ipo a 4>v$, (14) 

J R3 2 

4 = ~ [ d 3 x xl d>i (povf) A <h ( 15 ) 

4 = ~ J R3 d 3 x \ U 'i (pi + 2 Po + Po (4 0 - 2t/ 0 + 2r 2 (4’) 2 )) , (16) 

4 = f^d 3 xx l diU[ (pi + 2po + Po(*o-2t/ 0 + 2r 2 (4) 2 )). (17) 

In Table Q] we report results of the convergence tests for a sample system. In our 
implementation, numerical precision is controlled by the resolution of the grid, the maximum 
number L of Legendre polynomials used in the angular expansion of the solutions of the scalar 
and vector Poisson equations, and a value of the maximal difference between density distributions 
obtained in the last two consecutive iterations p to (In each iteration we compute the quantity 
Perr = maxjj \pf'j +l> — pfj\- Here index k numbers subsequent iterations; indices i and j refer 
to different grid nodes. The iteration procedure is stopped, when p err < p to i.) 

The results of the virial test do depend on the grid resolution and almost do not depend on 
the numbers of Legendre polynomials and the precision in convergence procedure. 


Table 1. Typical dependence of the results on the resolution of the numerical grid, the maximum 
number of the Legendre polynomials L, and the tolerance coefficient p to These results are 
obtained for a polytropic fluid with polytropic index 7 = 5/3 and the Keplerian rotation law 
Vq = wo/r 3 / 2 . The mass of the disk is = 0.327M c , the inner and outer radii are r m = 50i?g 

and r ou t = 500i?s respectively. 


Resolution 

L 

Ptol 

e 

e' 

e" 

100 

X 

100 

100 

10 -b 

2.53 x 10” b 

2.29 x 10' 4 

3.51 x 10" 4 

200 

X 

200 

100 

10 -6 

6.36 x 10 -6 

5.87 x 10~ 5 

8.87 x 10 -5 

400 

X 

400 

100 

10 -6 

1.56 x 10” 6 

1.93 x 10~ 5 

2.29 x 10" 5 

800 

X 

800 

100 

10 -6 

3.66 x 10" 7 

3.65 x 10~ 6 

4.62 x 10 -6 

1200 x 1200 

100 

10 -6 

1.44 x 10~ 7 

1.56 x 10” 6 

2.49 x 10 -6 

1600 x 1600 

100 

10 -6 

6.66 x 10 -8 

5.43 x 10“ 7 

1.39 x 10 -6 

400 

X 

400 

50 

10 -6 

1.56 x 10“ 6 

1.93 x 10~ 5 

2.27 x 10“ 5 

400 

X 

400 

75 

10 -6 

1.56 x 10~ 6 

1.93 x 10~ 5 

2.29 x 10" 5 

400 

X 

400 

100 

10 -6 

1.56 x 10 -6 

1.93 x 10~ 5 

2.29 x 10 -5 

400 

X 

400 

125 

10 -6 

1.56 x 10 -6 

1.93 x 10“ 5 

2.30 x 10 -5 

400 

X 

400 

150 

10 -6 

1.56 x 10“ 6 

1.93 x 10~ 5 

2.30 x 10 -5 

400 

X 

400 

100 

10 -5 

1.52 x 10~ 6 

1.94 x 10' 5 

2.29 x 10" 5 

400 

X 

400 

100 

10 -6 

1.56 x 10~ 6 

1.93 x 10~ 5 

2.29 x 10" 5 

400 

X 

400 

100 

10- 7 

1.59 x 10~ 6 

1.93 x 10- 5 

2.29 x 10 -5 

400 

X 

400 

100 

10 -8 

1.59 x 10 -6 

1.93 x 10~ 5 

2.29 x 10 -5 
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